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We calculate angular correlation function between ultra-high energy cosmic rays (UHECR) observed by 
Yakutsk and AGASA experiments, and most powerful BL Lacertae objects. We find significant correlations 
which correspond to the probability of statistical fluctuation less than 10~*, including penatly for selecting 
the subset of brightest BL Lacs. We conclude that some of BL Lacs are sources of the observed UHECR and 
present a list of most probable candidates. 
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Introduction. Identification of sources of ultra-high 
energy cosmic rays (UHECR) is extremely important. 
■Knowing production sites of UHECR will help to ex- 
plain the apparent absence of the Greisen-Zatsepin- 
Kuzmin (GZK) cutoff by selecting a particular class 
,of models. In the case of astrophysical origin it will 
give an invaluable information on physical conditions 
and mechanisms which may lead to acceleration of par- 
otides to energies of order 10^° eV. In the case of ex- 
tragalactic origin, it will provide a direct information 
about poorly known parameters which influence propa- 
gation of UHECR, such as extragalactic magnetic fields 
and universal radio background. 

There are observational reasons to believe that 
UHECR are produced by compact sources. It has been 
known for quite a while that the observed highest energy 
cosmic rays contain doublets and triplets of events com- 
ing from close directions |[ |^. Our recent analysis 
1^ based on the calculation of angular correlation func- 
tion shows that explanation of clusters by chance coin- 
cidence is highly improbable: the correlation function 
for Yakutsk events with energies E > 2.4 x 10^^ eV 
has an excess at 4° which would occur with probability 
2 X 10^^ for the uniform distribution, while the corre- 
lation function for AGASA events Q with energies 
E > 4.8 X 10^^ eV has an excess at 2.5° corresponding 
to chance probability 3 x 10^**. The combined proba- 
bility of the fluctuation in both sets is 4 x 10~^. So 
signiflcant autocorrelations should imply also large cor- 
relation of these events with their actual sources. It is a 
purpose of the present paper to identify these sources. 

The clustering of UHECR by itself imposes certain 
constraints on possible source candidates. With the ob- 
served fraction of events in clusters, the total number 



of sources can be estimated along the lines of Ref. Q to 
be of order several hundred. If the GZK cutoff is absent 
(or at energies below the cutoff), this estimate gives the 
number of sources in the entire Universe. Thus, to pro- 
duce observed clustering, the extragalactic sources have 
to be extremely rare as compared to ordinary galaxies. 
Taking 10'^ uniformly distributed sources for an esti- 
mate, the closest one is at z ~ 0.1. 

Various astrophysical candidates such as neutron 
stars, supernovae, gamma-ray bursts, colliding galax- 
ies, active galactic nuclei (AGN), lobes of radio-galaxies, 
dead quasars and others (for a review see Refs.|^ and 
references therein) have been proposed as sources of 
UHECR. Possible connection of highest-energy cosmic 
rays with these objects was considered in Refs. ^, [l0|. 
In this paper we study correlations of UHECR with 
BL Lacertae (BL Lac) objects which comprise a sub- 
class of AGN. Our motivations for selecting BL Lacs 
are as follows. If AGNs are sources, only those which 
have jets directed along the line of sight, or blazars, can 
correlate with observed UHECR events (regardless of 
the distance to a blazar in a world without GZK cut- 
off), since particles accelerated in a relativistic jet are 
strongly beamed. Blazars include BL Lacs and violently 
variable quasars with flat and highly polarized spectra. 
These spectral features give direct indication of seeing a 
relativistically beamed jet very close to the line of sight. 
BL Lacs is a subclass of blazars characterized, in ad- 
dition to the above spectral features which they share, 
by the (near) absence of emission lines in the spectra. 
This very important distinction indicates low density of 
ambient matter and radiation and, therefore, more fa- 
vorable conditions for acceleration to highest energies. 

The most recent catalog of AGNs and quasars con- 
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tains 306 confirmed BL Lacs [|Tl| . While this is the rich- 
est catalog we are aware of, it still may be incomplete. 
However, this is not crucial for establishing correlations 
between BL Lacs and UHECR events. Correlations of 
BL Lacs with UHECR were not studied before. We 
show that these correlations do exist and are statisti- 
cally significant. 

Method and results. Our method is based on calcu- 
lation of the angular correlation function and is similar 
to the one we have used in Ref. j|] . For each BL Lac, we 
divide the sphere into concentric rings (bins) with fixed 
angular size. We count the number of events falling into 
each bin and then sum over all BL Lacs, thus obtaining 
the numbers Ni (data counts). We repeat the same pro- 
cedure for a large number (typically 10^) of randomly 
generated sets of UHECR events. This gives the mean 
Monte-Carlo counts Nf^'^, the variance crf^^ and the 
probability p{5) to match or exceed the data count ob- 
served in the first bin. This probability is a function 
of the bin size S. Peaks of {Ni — Nf^'^)/ai or minima 
of p{5) with respect to 5 show angular scales at which 
correlations are most significant. 

The Monte-Carlo events are generated in the horizon 
reference frame with the geometrical acceptance 

dn oc cos Oz sin0 zd0z, 

where dz is the zenith angle. Coordinates of the events 
are then transformed into the equatorial frame assum- 
ing random arrival time. This transformation depends 
on the latitude of the experiment, so events simulat- 
ing different experiments are generated separately. The 
distribution of the generated Monte-Carlo events in dec- 
lination and right ascension reproduces well that of the 
experimental data. 

We have shown in Ref. ^ that autocorrelations are 
most significant for the two sets of UHECR events: 26 
Yakutsk events with energy E > 2.4 x 10^^ eV and 
39 AGASA events with energy E > 4.8 x 10^^ eV. If 
BL Lacs are sources of UHECR, their correlations with 
UHECR should be particularly large for these two sets. 
Assuming that energies of the events are not important 
for correlations at small angles, we combine them to- 
gether in one set of 65 events. 

Since acceleration of particles to energies of order 
10^° eV typically requires extreme values of parame- 
ters, probably not all BL Lacs emit UHECR of required 
energy. We assume that this ability is correlated with 
optical and radio emissions, and select the most pow- 
erful BL Lacs by imposing cuts on redshift, apparent 
magnitude and 6 cm radio flux. For more than a half 
of BL Lacs the redshift is not known. It is generally 
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FIG. 1. The sky map (in Galactic coordinates) with 65 
UHECR events (circles) and BL Lacertae objects with 
cuts (|l|). 

Table 1. Names and coordinates (Galactic longitude, 
latitude and redshift) of BL Lacs plotted in Fig. ^ which 
fall within 3° from some UHECR event (their energies 
are listed in the last column). 



Name 


1 


b 


z 


E/IO^^ eV 


lES 0806+524 


166.25 


32.91 


0.138 


3.4; 2.8; 2.5 


RX J10586+5628 


149.59 


54.42 


0.144 


7.76; 5.35 


2EG J0432+2910 


170.52 


-12.6 




5.47; 4.89 


OT 465 


74.22 


31.4 




4.88 


TEX 1428-1-370 


63.95 


66.92 


0.564 


4.97 



expected that these BL Lacs are at z > 0.2. We include 
them in the set. The cuts 

z > 0.1 or unknown; mag < 18; Fq > 0.17 Jy (1) 

leave 22 BL Lacs which are shown in Fig. Q together with 
65 cosmic rays from the combined set. The dependence 
on cuts is discussed below. 

As one can see from Fig. |l|, two of 22 BL Lacs coin- 
cide with the two triplets of UHECR events, one coin- 
cides with a doublet and two BL Lacs lie close to sin- 
gle events. This is reflected in the correlation function, 
which is plotted in Fig. |^ for the bin size 2.5°. It has 
8 events in the first bin while 1.25 is expected for the 
uniform distribution. The probability of such an excess 
is 2 X 10^^. BL Lacs and UHECR events which con- 
tribute to this correlation are listed in Table 1. Note 
that at large angles the correlation function fluctuates 
around zero, which shows that the acceptance in the 
Monte-Carlo simulation is chosen correctly. 

The probability p{S) as a function of the angular 
separation S is shown in Fig. ^. It has a minimum at 
2.5°. For comparison, smooth curve shows the behavior 
expected when 9 events out of 65 come from BL Lacs 
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FIG. 2. The angular correlation function between the 
combined set of UHECR and BL Lac set (|l|). 



10-1 



■.^ 10-2 r 



10-3 



10-" 



10-5 




1 2 3 4 5 6 7 8 9 10 1 1 12 13 14 15 16 
bin size/l° 

FIG. 3. The dependence of the probability p{S) on the 
bin size 5 for the combined set of UHECR and BL Lac 
set (§). 



(assuming that accuracy of angle determination is 1.8° 
and distribution of errors is Gaussian ). 

The smah angular size of the peak in the correla- 
tion function, compatible with the experimental angular 
resolution, suggests that UHECR events responsible for 
these correlations are produced by neutral primary par- 
ticles. Indeed, if the primaries were charged they would 
have been deflected in the Galactic magnetic field by 
3° — 7° depending on arrival direction, particle energy 
and the model of the magnetic field, and correlations at 
2.5° would be destroyed. 

Discussion. We have seen that 22 bright BL Lacs and 
65 cosmic rays from the combined set are strongly cor- 
related: the probability to find 8 or more out of 65 ran- 
domly generated cosmic rays within 2.5° of any of the 
BL Lacs is 2 X 10~^. Should one conclude that BL Lacs 
are sources of UHECR, or the above correlation may be 



an artifact of our selection procedure? Let us discuss 
possible loopholes. 

First potential source of problem is incompleteness 
of the BL Lac catalog and non-uniform coverage of the 
sky. Indeed, 22 BL Lacs selected by cuts (1) almost all 
lie in the Northen hemisphere due to observational bias. 
However, it is easy to understand that, unlike for many 
other astrophysical problems, for establishing the fact 
of correlations with UHECR the incompleteness of BL 
Lac catalog is not essential. The method we use works 
for any set of potential sources regardless of their dis- 
tribution over the sky (including such extreme cases as 
just one source, or a compact group of several sources). 
This is guaranteed by using the same set of sources with 
real data and with each Monte-Carlo configuration. 

Second potential problem is related to the fact that 
there exist strong autocorrelations in the UHECR set, 
while Monte-Carlo events are not correlated. One may 
wonder if the observed correlation with BL Lacs is (par- 
tially) due to autocorrelations of UHECR. To see that 
this effect is negligible in our case, we performed test 
Monte-Carlo simulations with configurations containing 
the same number of doublets and triplets as the real 
data, and random in other respects. We found practi- 
cally no difference between the two methods. 

Finally, there is an issue of cuts and related issue of 
selection of catalogs. One may worry that by adjusting 
several cuts and searching in several catalogs the prob- 
ability as small as Pmin ~ 10~^ can be found with any 
set of astrophysical objects, even with those which have 
nothing to do with UHECR. So, the question is how 
easily the low values of Pmin can be obtained within the 
adopted procedure of cuts. This question can be stud- 
ied quantitatively by assigning a proper penalty for each 
try in such a way that resulting probability gives true 
measure for the correlations in question to be a statis- 
tical fluctuation. For the case at hand we have found 
that when proper penalties are assigned, the resulting 
probability is larger than Pmin by about an order of mag- 
nitude. In other words, one would have to try thousands 
of catalogs to find correlation as significant as we have 
found for BL Lacs. We present the procedure of penalty 
calculation and resulting significance of correlations be- 
low. 

In fact, we did not search for correlations with other 
catalogs of astrophysical objects. Thus, no penalty is 
associated with that. Similarly, we did not adjust the 
set of cosmic rays (as explained before, it was selected 
in Rcf. 1^ on the basis of most significant autocorre- 
lations). But we do adjust cuts in the BL Lac catalog. 
Therefore, we should assign a penalty factor to this ad- 
justment. 
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It is clear that some cuts have to be made because 
65 events may have at most 65 sources among 306 BL 
Lacs in the catalog (probably much less). In our cal- 
culations we imposed cuts on redshift, magnitude and 
6 cm radio-flux. The cut on redshift is motivated by the 
expected total number of sources; we did not adjust this 
cut to minimize the probability. Cuts on magnitude and 
radio-flux were adjusted. Corresponding penalty can be 
calculated in the following way (cf. Ref d). A rando m 
set of cosmic rays should be generated and treated as 
real data, i.e. minimum probability Pmin is searched for 
by adjusting the cuts in the BL Lac catalog in exactly 
the same way as it was done for the real data. This 
should be repeated many times, giving different Pmin 
each time. The number of occurrences of a given value 
of Pmin is then calculated as a function of Pmin- This 
gives the probability (we call it Pcor) that the adjust- 
ment of the cuts in BL Lac catalog produces p < Pmin 
with a random set of cosmic rays. The probability Pcor 
is a correct measure of the significance of correlations. 
We define Pcor/Pmin > 1 as the penalty factor. 

We calculated pcor with 10^ random sets of cos- 
mic rays. We have found that the penalty grows at 
small Pmin and approaches a constant value in the limit 
Pmin ^ (for this reason it is more convenient to define 
the penalty factor than to work in terms of Pcor)- For 
the real set of UHECR Pmin = 4 x 10^^ and is reached 
with the cuts 

z > 0.1 or unknown; mag < 16; Fq > 0.17 Jy (2) 

They leave 5 BL Lacs two of which coincide with 
triplets. (In the previous section different cuts are pre- 
sented because, with similar significance, they include 
more potential sources.) This probability should be 
multiplied by the penalty factor. We found that the 
penalty factor is ~ 15 at Pmin — 10"^ This gives 

Pcor = 6 X 10~^, which is the probability that the cor- 
relation we have found is a statistical fluctuation. 

Conclusions. The significant correlations between 
UHECR and BL Lacs imply that at least some of BL 
Lacs are sources of UHECR. Most probable candidates 
can be seen in Fig. ^ and are listed in Table 1. Two 
BL Lacs, lES 0806-1-524 and RX J10586-I-5628, coin- 
cide with triplets of UHECR events (in the second case 
the third event of a triplet is at 4.5° and is not listed 
in the table). Both of them are at the distance of 
~ 600 Mpc from the Earth. The next-probable can- 
didate 2EG J0432-h2910 has unknown redshift. 

The correlations at small angles arc difficult to ex- 
plain by charged primary particles. Within the Stan- 
dard Model the only two neutral candidates are photon 



and neutrino. Photon attenuation length at i? < 10^° 
eV is much smaller (see e.g. [^) than the distance to 
even the closest BL Lac. However, photons can not be 
ruled out yet if one assumes sources at d ^ 600 Mpc 
and "extreme" astrophysical conditions: primary parti- 
cles accelerated to _E > 10^^ eV with "hard" spectrum 
~ E~" and a < 2, and extragalactic magnetic fields 
B < 10"" G |l|]. Neutrino models || require simi- 
lar assumptions except that constraints on the magnetic 
filed are relaxed for "pure" neutrino sources and there is 
no constraint on the distance to the sources. However, if 
"pure" neutrino sources cannot be arranged, the model 
effectively becomes "photonic" ||l^ . If astrophysical dif- 
ficulties can be overcome, these models will be appealing 
candidates for the solution of the UHECR puzzle. Alter- 
natively, one may resort to a new physics, e.g., violation 
of the Lorentz invariance p^ . 

Independent cross-checks are necessary to determine 
whether particular objects are sources of UHECR. One 
of these cross-checks could be coincidence of arrival time 
of events contributing to small angle correlations with 
periods of activity of candidate BL Lacs. Dedicated 
monitoring of these BL Lac may be suggested. It is also 
important to analyze possible specific properties of air 
showers initiated by these events. 
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